

#include <vector.h>
#include <math.h>

#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <stdlib.h>
#include "molCluster.h"

using namespace std;

molCluster::molCluster(){

	nMolecules = 0;
	nAtoms = 0;
	nUniqueBonds = 0;
}

molCluster::molCluster( const molCluster& rhs ) {
	int i;

        nMolecules = rhs.nMolecules;
        nAtoms = rhs.nAtoms;

	atoms.resize( nAtoms );
	for( i=0; i<nAtoms; i++ ) {
        	atoms[i] = rhs.atoms[i];
	}


        properties = rhs.properties;

	for( i=0; i<SIGLENGTH; i++ ) {
		signature[ i ] = rhs.signature[ i ]; 
	}
	
        molConnectList = rhs.molConnectList;

        molLabelList = rhs.molLabelList;

	evaluateFunc = rhs.evaluateFunc;

	nUniqueBonds = rhs.nUniqueBonds;

	molBondList = rhs.molBondList;
}

molCluster::molCluster( string& molSettingFileName, vector<double>(*anEvaluator)(vector<atom>& a, int optimize) ) {
	
	int i, j, nRow;

	nAtoms = 0;

	ifstream settingFile( molSettingFileName.c_str() );

	string dummy;

	settingFile >> dummy >>  nRow >> dummy;
	vector<int> nCol(nRow);
	for(i=0; i<nRow; i++) settingFile >> nCol[i];


	molConnectList.resize( nRow );
	settingFile >> dummy;
	for( i=0; i<nRow; i++ ) {
		molConnectList[ i ].resize( nCol[i] );
		for( j=0; j<nCol[i]; j++ ) {
			nAtoms++;
			settingFile >> molConnectList[ i ][ j ];
		}
	}	


	settingFile >> dummy;
	molLabelList.resize( nRow );
	for( i=0; i<nRow; i++ ) {
		molLabelList[ i ].resize( nCol[i] );
		for( j=0; j<nCol[i]; j++ ) {
			settingFile >> molLabelList[ i  ][ j ];
		}
	}
	
	settingFile >> dummy >> nUniqueBonds >> dummy;
	molBondList.resize(nUniqueBonds);
	for( i=0; i<nUniqueBonds; i++ ) {
		settingFile >> molBondList[i].label1 >> molBondList[i].label2 >> molBondList[i].maxDistance;
	}

	settingFile.close();

	nMolecules = nRow;
	atoms.resize( nAtoms );
	properties.resize( NPROPERTIES );

	evaluateFunc = anEvaluator;
}

molCluster::~molCluster( ){ }	

//gets the number of atoms
int molCluster::size( ) {
	return atoms.size();
}
			
	
//get the centroid
void molCluster::getCentroid( atom &centroid ) {

	centroid.x = 0;
	centroid.y = 0;
	centroid.z = 0;
		
	//compute the centroid
	for (int j=0; j<3; j++) {
		for (int i=0; i<atoms.size(); i++) {
			centroid.setAt(j, centroid[j] + atoms[i][j]);
		}
		centroid.setAt(j, centroid[j]/(double) atoms.size());
	}	
}
	
//set to centroid;
void molCluster::setToCentroid(void) {
	atom centroid;
	getCentroid(centroid);
		
	//set to the centroid;
	for (int i=0; i<atoms.size(); i++)
		atoms[i] -= centroid;

	return;
}
	
/*
   Rotate a point p by angle theta around an arbitrary axis r
   Return the rotated point.
   Positive angles are anticlockwise looking down the axis
   towards the origin.
   Assume right hand coordinate system.
*/
void molCluster::ArbitraryRotate(atom p, atom &q, double theta, atom r) {
	q = p;
	normalise(r);
	double q0 = cos(theta/2);
	double q1 = sin(theta/2)*r.x;
	double q2 = sin(theta/2)*r.y;
	double q3 = sin(theta/2)*r.z;

	q.x = p.x*(q0*q0 + q1*q1 - q2*q2 - q3*q3) + p.y*(2*(q1*q2 - q0*q3)) + p.z*(2*(q1*q3 + q0*q2));
	q.y = p.x*(2*(q2*q1 + q0*q3)) + p.y*(q0*q0 - q1*q1 + q2*q2 - q3*q3) + p.z*(2*(q2*q3 - q0*q1));
	q.z = p.x*(2*(q3*q1 - q0*q2)) + p.y*(2*(q3*q2 + q0*q1)) + p.z*(q0*q0 - q1*q1 - q2*q2 + q3*q3);
}		

void molCluster::normalise(atom &p) {
	double s;
	s = sqrt( (p.x * p.x) + (p.y * p.y) + (p.z * p.z));
	p.x = p.x/s;
	p.y = p.y/s;
	p.z = p.z/s;
}

//rotates a cluster about the centroid by theta (in radians)
//WARNING: Assumes that the centroid is centered around (0,0,0)
void molCluster::rotateAboutCentroid(double theta) {
	//first find an arbitrary point in the space, R
	//do this by randomly translating the centroid
	atom R;
	R.x = rand();
	R.y = rand();
	R.z = rand();
		
	for (unsigned int i=0; i<atoms.size(); i++) {
		atom temp;
		ArbitraryRotate(atoms[i], temp, theta, R);
		atoms[i] = temp;
	}
}	
	
	
//save to ofstream
void molCluster::printToFileStream( ofstream& out ) {

	//output size
	out << scientific << setprecision( PRECISION );
	out << size() << endl;
	
	//output properties
	out << properties.size() << " ";
	for (int i=0; i<properties.size(); i++) {
		out << properties[i] << " ";
	}
	out << endl;
		
	for (int i=0; i<size(); i++) {
		atoms[ i ].printToFileStream( out );
	}

}
	
	
//assignment operator
molCluster& molCluster::operator=( const molCluster &rhs) {
	atoms.clear();
	atoms.resize(rhs.nAtoms);
	for (int i=0; i<rhs.nAtoms; i++) {
		atoms[i] = rhs.atoms[i];
	}
		
	properties = rhs.properties;
		
	for (int i=0; i< SIGLENGTH; i++) {
		signature[i] = rhs.signature[i];
	}
		
	molConnectList = rhs.molConnectList;
	molLabelList = rhs.molLabelList;
	nAtoms = rhs.nAtoms;
	nMolecules = rhs.nMolecules;
		
	evaluateFunc = rhs.evaluateFunc;

	nUniqueBonds = rhs.nUniqueBonds;
	molBondList = rhs.molBondList;

	return *this;
}
	
//get a particular atom
atom& molCluster::operator[](int index) {
	if (index >= size() || index < 0) {
		cout << "ERROR! trying to access non-existent entry." << endl;
	}
	return (atoms[index]);
}
	
	
//Ultrafast Shape Recognition Signature computations

/*	gets statistics from a point/atom a
 *	computes the mean, standard deviation and assymetry of distances from a
 *	i.e. the first, second and 3rd standardized moments.
 *	returns false if the cluster is smaller than 2 atoms since no statistics can be derived
 */
bool molCluster::getStatistics( atom & a , double &mean, double &sd, double &assymetry ) {
	int i;
	int size = atoms.size();
		
	if (size <= 1) return false;
		
	//compute mean
	mean = 0;
	for (i=0; i<size; i++) {
		mean += a.getDistanceTo(atoms[i]);
	}
	mean /= size;
		
	//compute standard deviation
	sd = 0;
	for (i=0; i<size; i++) {
		sd += pow(a.getDistanceTo(atoms[i]) - mean, 2.0);
	}
		
	sd = sqrt(sd/(double) (size - 1.0));
		
	//compute assymetry/skewness
		
	double thirdmoment = 0;
	for (i=0; i<size; i++) {
		thirdmoment += pow((a.getDistanceTo(atoms[i]) - mean), 3.0);
	}
	if (size > 2) {
		thirdmoment /= (double) ((size - 1)*(size-2));
		thirdmoment = thirdmoment*size;
	}
	else {
		thirdmoment /= (double) (size - 1);	
	}
	assymetry = thirdmoment/pow(sd, 3.0);
		
	return true;
}

//calculates the USR signature associated with this cluster
bool molCluster::calcSignature(void) {

	//create the reference atoms
	atom centroid;
	getCentroid(centroid);
	atom catom;
	getClosestAtomTo(centroid, catom);
	atom fatom;
	getFurthestAtomFrom(centroid, fatom);
	atom ffatom;
	getFurthestAtomFrom(fatom, ffatom);

	//set signatures to zero
	//resetSignature();

	//calculate mean, std, assymetry of distances from centroid
	int errorCode = true;
	errorCode = errorCode && (getStatistics(centroid, signature[0], signature[1], signature[2]));
	errorCode = errorCode && (getStatistics(catom, signature[3], signature[4], signature[5]));
	errorCode = errorCode && (getStatistics(fatom, signature[6], signature[7], signature[8]));
	errorCode = errorCode && (getStatistics(ffatom, signature[9], signature[10], signature[11]));
		
	return errorCode;
}

//gets the USR signature and puts it into sig
//all items currently in sig are removed
//returns false if size of cluster is 1 or less.
bool molCluster::getSignature(vector <double> &sig) {
	sig.clear();
	int i;
	for (i=0; i<SIGLENGTH; i++) {
		sig.push_back(signature[i]);
	}
	return (atoms.size() > 1);
}

//resets the internal signature
void molCluster::resetSignature( ) {
	int i;
	for (i=0; i<SIGLENGTH; i++) signature[i] = 0;
	return;
}

//calculates the similarity to rhs
//uses the USR metric
double molCluster::getSimilarityTo( molCluster &rhs ) {
	double S = 0;
	int i;
		
	//compute manhattan distance between the two signatures
	for (i=0; i<SIGLENGTH; i++) {
		S += fabs(signature[i] - rhs.signature[i]);
	}
	//get the inverse and scaled metric
	S /= 12.0;
	S = 1.0/(1+S);
	return S;	
}

//returns the euclidean distance to rhs's signature
double molCluster::getSigEucDistanceTo(molCluster &rhs) {
	double S = 0;
	int i;
		
	//compute manhattan distance between the two signatures
	for (i=0; i<SIGLENGTH; i++) {
		S += pow(signature[i] - rhs.signature[i], 2.0);
	}
	S = sqrt(S);
		
	return S;	
}	

//for external calls to evaluate any cluster of atoms and not only my own
int molCluster::getFurthestAtomFrom( atom &c, atom &fatom) {
	int i;
	int size = atoms.size();
		
	//get the furthest atom
	int furthestAtom = 0;
	double furthestDist = atoms[0].getDistanceTo(c) ;
	double currentDist = 0;
	for (i=1; i<size; i++) {
		currentDist = atoms[i].getDistanceTo(c);
		if (currentDist > furthestDist) {
			furthestAtom = i;
			furthestDist = currentDist;
		}
	}
		
	fatom = atoms[furthestAtom];
	return furthestAtom;
}

int molCluster::getClosestAtomTo( atom &c, atom &catom) {
	int i;
	int size = atoms.size();
		
		
	//get the closest atom
	int closestAtom = 0;
	double closestDist = atoms[0].getDistanceTo(c) ;
	double currentDist = 0;
	for (i=1; i<size; i++) {
		currentDist = atoms[i].getDistanceTo(c);
		if (currentDist < closestDist) {
			closestAtom = i;
			closestDist = currentDist;
		}
	}
		
	catom = atoms[closestAtom];
	return closestAtom;
		
}

void molCluster::mutate( int strength ) {


        double test = mt();

	vector<atom> tempAtoms = atoms;

	int nAttempts = 0;

	while(nAttempts < MAX_MUTATION_ATTEMPTS) {
        	if (test < 0.8) {

        		//in general, we want a small mutation
        		for (int i=0; i< strength; i++) {
                		//randomly pick a molecule
                        	int pickedMol = ((int) (mt()*10000))%( nMolecules );

                        	//perform a translation
                        	molTranslation( pickedMol, 1.0 ); //translate by a max of 2 angstroms

                        	//perform a rotation
                        	molRotation( pickedMol );
                 	}
         	} else if (test < 0.9) {


         		//and sometimes, we want a mix
                	//but sometimes we want a more destructive mutation
                	for (int i=0; i< strength; i++) {
                		//randomly pick a molecule
                        	int pickedMol = ((int) (mt()*10000))%( nMolecules );

                        	test = mt();
                    	

				if (test < 0.5) {
                 
			       		//relocate a random molecule to the surface
                                	relocateMol( pickedMol );
                        	} else {
				
                        		//perform a translation
                                	molTranslation( pickedMol, 1.0 ); //translate by a max of 2 angstroms

                                	//perform a rotation
                                	molRotation( pickedMol );
                        	}
                	}

          	} else {


          		//but sometimes we want a more destructive mutation
                	for (int i=0; i< strength; i++) {
                		//randomly pick a molecule
                        	int pickedMol = ((int) (mt()*10000))%( nMolecules );

                        	//relocate a random molecule to the surface
                        	relocateMol( pickedMol );
                 	}
         	}

		
		if( isValid( ) )  break; 
		else atoms = tempAtoms;

		nAttempts++;
	}

        //compute USR signature
        calcSignature();

}



//mutation types
//just translate a molecule in c by a distance (limited by maxAxisDistance) along each axis
void molCluster::molTranslation( int pickedMol, double maxAxisDistance ) {

        //create a translation vector
        atom translation;
        for ( int axis = 0; axis < 3 ; axis++ ) {
        	translation.setAt( axis, mt()*maxAxisDistance );
        }

	// translate the atoms of pickedMol
	int nAtomsInPickedMol = molConnectList[pickedMol].size();
	for( int i=0; i<nAtomsInPickedMol; i++ ) {
		atoms[ molConnectList[ pickedMol ][ i ]  ] += translation;
	}

}

//rotate a molecule around its center of mass
void molCluster::molRotation( int pickedMol ) {

       /*
              rotating a bound hydrogen atom randomly is equivalent to picking a point on the sphere (with the same
              distance as the original hydrogen atom) and creating a new hydrogen atom there (and destroying the 
              original).
        */

	int nAtomsInPickedMol = molConnectList[pickedMol].size();
	for( int i=1; i<nAtomsInPickedMol; i++ ) {
		double distance = atoms[  molConnectList[pickedMol][i]  ].getDistanceTo(  atoms[ molConnectList[pickedMol][0]  ]   );
               	double u = mt()*2.0 -1.0;
               	double theta = mt()*2*M_PI;

		atom translation;
        	translation.x = distance*sqrt(1.0 - u*u)*cos(theta);
                translation.y = distance*sqrt(1.0 - u*u)*sin(theta);
                translation.z = distance*u;
	
		translation += atoms[ molConnectList[pickedMol][0]  ];
		translation.label = molLabelList[ pickedMol ][ i ];

		atoms[  molConnectList[pickedMol][i]  ] = translation;
	}

}

//relocate molecule 
void molCluster::relocateMol( int pickedMol ) {

	/*
                Moving the a molecule (m) to another spot on the surface of the cluster in 4 easy steps:
                1. Find d, the distance from f (the outermost mol) to the centroid
                2. Pick p, a point on a sphere with distance d from the centroid
                3. Calculate v, the translation vector that will move f to p. 
                4. Translate each of f's atoms.
         */
         atom centroid;
         getCentroid(centroid);

         //find maximum distance from centroid
         double maxDistance = -1;
         for (int i=0; i<size(); i++) {
         	if ( atoms[i].getDistanceTo(centroid) > maxDistance ) {
                	maxDistance = atoms[i].getDistanceTo(centroid);
                }
         }

	//
	//	?? is this value of 2.0 is general for all ??
	//
        maxDistance += 2.0; //we add 2.0 Angstroms to be sure of no overlap
        double u = mt()*2.0 -1.0;
        double theta = mt()*2*M_PI;

        atom translation;
        translation.x = maxDistance*sqrt(1.0 - u*u)*cos(theta);
        translation.y = maxDistance*sqrt(1.0 - u*u)*sin(theta);
        translation.z = maxDistance*u;


	int nAtomsInPickedMol = molConnectList[pickedMol].size();
        for( int i=0; i<nAtomsInPickedMol; i++ ) {
		atoms[  molConnectList[pickedMol][i]  ] += translation;
	}

}



void molCluster::crossover( molCluster mum, molCluster& child ) {

	int i, j;
	
	molCluster dad = (*this);
	

	vector<atom> latestAtoms(nAtoms);
	vector<atom> tempAtoms = atoms;

	molCluster tempCluster = (*this);

	int index;
	for(i=0; i<nMolecules; i++) {
		
		latestAtoms = tempAtoms;
		
		//try use dad
		if(mt()<0.5) {
			for(j=0; j<molConnectList[i].size(); j++) {
				index = molConnectList[i][j];		
				tempAtoms[ index  ] = dad.atoms[ index ];
			}
			// if failed, then use mum
			tempCluster.atoms = tempAtoms;
			if( !tempCluster.isValid( )) {
				for(j=0; j<molConnectList[i].size(); j++) {
                        		index = molConnectList[i][j];
                                	tempAtoms[ index  ] = mum.atoms[ index ];
				}
			}

			tempCluster.atoms = tempAtoms;
			if(!tempCluster.isValid( )) {
				tempAtoms = latestAtoms;
			}	
				
			
		}
		else {	// use mum


			for(j=0; j<molConnectList[i].size(); j++) {
                        
				index = molConnectList[i][j];
                                tempAtoms[ index  ] = mum.atoms[ index ];

			}
			// if failed, then use dad
			tempCluster.atoms = tempAtoms;
                        if( !tempCluster.isValid( )) {
                                for(j=0; j<molConnectList[i].size(); j++) {
                                
					index = molConnectList[i][j];
                                	tempAtoms[ index  ] = dad.atoms[ index ];

				}
                        }

			tempCluster.atoms = tempAtoms;
                        if(!tempCluster.isValid( )) {
                                tempAtoms = latestAtoms;
                        }


		}

		
	}	

	child = dad;
	child.atoms = tempAtoms;
}


void molCluster::evaluate( int optimize ) {


	vector<double> output(NPROPERTIES);

	output = evaluateFunc( atoms, optimize );

	for(int i=0; i<NPROPERTIES; i++) {
		properties[ i ] = output[ i ];
	}

}

bool molCluster::molCluster::isValid() {

        int i, j;

        double distance;

        vector< vector<int> > d(nAtoms);
        for(i=0; i<nAtoms; i++) d[i].resize(nAtoms);

        for(i=0; i<nAtoms; i++) {
                for(j=i; j<nAtoms; j++) {
                        distance = atoms[i].getDistanceTo(atoms[j]);
                        if( checkBondConnect( atoms[i].label, atoms[j].label, distance )  ) {
                                d[i][j]=1;
                                d[j][i]=1;
                        }
                }
        }


        return bConnectivity( d, nAtoms );

}

bool molCluster::bConnectivity(vector<vector<int> > d,int n){
        int i, j, l;
        for( i=0; i<n; i++){
                for( j=0; j<n; j++ )
                        if( (i!=j) && (d[i][j]==1) ) {//if connected,jump to jth row and flip
                                for( l=0; l<n; l++) if(d[i][l]==1) d[j][l]=1;
                        }
        }
        for(i=0;i<n;i++)
                for(j=0;j<n;j++) if(d[i][j]==0) return false;

        return true;
}

bool molCluster::checkBondConnect( char a, char b, double distance  ) {

	for(int i=0; i<nUniqueBonds; i++) {
		bond& c = molBondList[i];
		if( c.label1==a && c.label2==b && distance<c.maxDistance ) return true;
		if( c.label1==b && c.label2==a && distance<c.maxDistance ) return true;
	}

        return false;
}

ostream& operator<<( ostream& out, molCluster& aMolCluster ) {

	//output size
        out << scientific << setprecision( PRECISION );
	out << aMolCluster.size() << endl;

        //output properties
        out << aMolCluster.properties.size() << " ";
        for (int i=0; i<aMolCluster.properties.size(); i++) {
                out << aMolCluster.properties[i] << " ";
        }
        out << endl;

        for (int i=0; i<aMolCluster.size(); i++) {
                out << aMolCluster.atoms[ i ];
        }
	return out;
}


